home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / prgtools / gnustuff / minix / libsrc~1.z / libsrc~1 / dflonum.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-12-28  |  14.6 KB  |  647 lines
  1. /* A working Gnulib68k, thanks to Scott McCauley for the origonal
  2.    version, and John Dunning (jrd@STONY-BROOK.SCRC.Symbolics.COM)
  3.    who got this working.
  4.    
  5.    Please not that only the double float. stuff is picked up from
  6.    here, the other long-int and single float stuff come from
  7.    fixnum.s and sflonum.s (see flonum.h for the appr. #defs).
  8.    ++jrb
  9.    */
  10.  
  11. /* Subroutines needed by GCC output code for the 68000/20. 
  12.    
  13.    Compile using -O flag with gcc. 
  14.    Use the -m68000 flag if you have a 68000
  15.    
  16.    This package includes 32 bit integer and 64-bit floating
  17.    point routines.
  18.    
  19.    WARNING: alpha-test version (July 1988) -- probably full of bugs.
  20.    If you find a bug, send bugs reports to jsm@phoenix.princeton.edu,
  21.    or
  22.    
  23.    Scott McCauley
  24.    PPPL P. O. Box 451
  25.    Princeton NJ 08543
  26.    
  27.    Known bugs/features:
  28.    
  29.    1) Depending on the version of GCC, this may produce a lot of
  30.    warnings about conversion clashes, etc. It should still compile
  31.    as is. Also, there appears to be a bug in the register-allocation
  32.    parts of gcc that makes the compiler sometimes core dump
  33.    or run into an infinite loop. This version works -- if you
  34.    decide to change routines these compiler bugs can bite very hard....
  35.    
  36.    2) all single-precision gets converted to double precision in the
  37.    math routines (in accordance with K&R), so doing things in
  38.    double precision is faster than single precision....
  39.  
  40. (not any more -- jrd )
  41.    
  42.    3) double precision floating point division may not be accurate to
  43.    the last four or so bits. Most other routines round to the
  44.    lsb.
  45.  
  46. (not any more -- jrd )
  47.    
  48.    4) no support of NaN and Inf
  49.    
  50.    5) beware of undefined operations: i.e. a float to integer conversion
  51.    when the float is larger than MAXIINT.
  52.    
  53.    */
  54.  
  55. #include "flonum.h"
  56.  
  57. #ifdef __GCC_OLD__
  58. #define umultl _umulsi
  59. #define multl _mulsi3
  60. #define udivl _udivsi3 
  61. #define divl _divsi3
  62. #define ddiv _divdf3
  63. #define qmult _qmult
  64. #define dmult _muldf3
  65. #define dneg _negdf2
  66. #define dadd _adddf3
  67. #define dcmp _cmpdf2
  68. #define dtoul _fixunsdfsi
  69. #define dtol _fixdfsi
  70. #define ltod _floatsidf
  71. #else
  72. #define umultl __umulsi
  73. #define multl __mulsi3
  74. #define udivl __udivsi3 
  75. #define divl __divsi3
  76. #define ddiv __divdf3
  77. #define qmult __qmult
  78. #define dmult __muldf3
  79. #define dneg __negdf2
  80. #define dadd __adddf3
  81. #define dcmp __cmpdf2
  82. #define dtoul __fixunsdfsi
  83. #define dtol __fixdfsi
  84. #define ltod __floatsidf
  85. #endif
  86.  
  87. #ifdef L_divdf3
  88.  
  89. /* new double-float divide routine, by jrd */
  90. /* thanks jrd !! */
  91.  
  92. #ifdef __GCC_OLD__
  93. double _divdf3(num, denom)
  94. #else
  95. double __divdf3(num, denom)
  96. #endif
  97. double num, denom;
  98. {
  99.     double local_num = num;
  100.     double local_denom = denom;
  101.     struct bitdouble * num_bdp = (struct bitdouble *)(&local_num);
  102.     struct bitdouble * den_bdp = (struct bitdouble *)(&local_denom);
  103.     short num_exp = num_bdp->exp,
  104.           den_exp = den_bdp->exp;
  105.     short result_sign = 0;
  106.     /*  register */ unsigned long num_m1, num_m2, num_m3, num_m4;
  107.     register unsigned long den_m1, den_m2, den_m3, den_m4;
  108.     unsigned long result_mant[3];
  109.     unsigned long result_mask;
  110.     short result_idx;
  111.     
  112.     if ((num_exp == 0) || (den_exp == 0)) /* zzz should really cause trap */
  113.       return(0.0);
  114.     
  115.     /* deal with signs */
  116.     result_sign = result_sign + num_bdp->sign - den_bdp->sign;
  117.     
  118.     /* unpack the numbers */
  119.     num_m1 = num_bdp->mant1 | 0x00100000;        /* hidden bit */
  120.     num_m2 = num_bdp->mant2;
  121.     num_m4 = num_m3 = 0;
  122.     den_m1 = den_bdp->mant1 | 0x00100000;        /* hidden bit */
  123.     den_m2 = den_bdp->mant2;
  124.     den_m4 = den_m3 = 0;
  125.     
  126. #if 0
  127.     /* buy a little extra accuracy by shifting num and denum up 10 bits */
  128.     for (result_idx /* loop counter */ = 0 ; result_idx < 10 ; result_idx++)
  129.     {
  130.     ASL3(num_m1, num_m2, num_m3);
  131.     ASL3(den_m1, den_m2, den_m3);
  132.     }
  133. #endif
  134.     
  135.     /* hot wire result mask and index */
  136.     result_mask = 0x00100000;            /* start at top mantissa bit */
  137.     result_idx = 0;                /* of first word */
  138.     result_mant[0] = 0;
  139.     result_mant[1] = 0;
  140.     result_mant[2] = 0;
  141.     
  142.     /* if denom is greater than num here, shift denom right one and dec num expt */
  143.     if (den_m1 < num_m1)
  144.       goto kludge1;                /* C is assembly language,
  145.                            remember? */
  146.     if (den_m1 > num_m1)
  147.       goto kludge0;
  148.     if (den_m2 <= num_m2)                /* first word eq, try 2nd */
  149.       goto kludge1;
  150.     
  151.   kludge0:
  152.     
  153.     num_exp--;
  154.     ASR4(den_m1, den_m2, den_m3, den_m4);
  155.     
  156.   kludge1:
  157.     
  158.     for ( ; !((result_idx == 2) && (result_mask & 0x40000000)) ; )
  159.     {
  160.     /* if num >= den, subtract den from num and set bit in result */
  161.     if (num_m1 > den_m1) goto kludge2;
  162.     if (num_m1 < den_m1) goto kludge3;
  163.     if (num_m2 > den_m2) goto kludge2;
  164.     if (num_m2 < den_m2) goto kludge3;
  165.     if (num_m3 > den_m3) goto kludge2;
  166.     if (num_m3 < den_m3) goto kludge3;
  167.     if (num_m4 < den_m4) goto kludge3;
  168.     
  169.       kludge2:
  170.     result_mant[result_idx] |= result_mask;
  171.     SUB4(den_m1, den_m2, den_m3, den_m4, num_m1, num_m2, num_m3, num_m4);
  172.       kludge3:
  173.     ASR4(den_m1, den_m2, den_m3, den_m4);
  174.     result_mask >>= 1;
  175.     if (result_mask == 0)        /* next word? */
  176.     {
  177.         result_mask = 0x80000000;
  178.         result_idx++;
  179.     }
  180.     }
  181.  
  182. /* stick in last half-bit if necessary */
  183.   if (result_mant[2])
  184.     {
  185.     den_m1 = 0;        /* handy register */
  186.     den_m2 = 1;
  187. /* busted
  188.     ADD2(den_m1, den_m2, result_mant[0], result_mant[1]);
  189. */
  190.     result_mant[1]++;
  191.     if (result_mant[1] == 0)
  192.         result_mant[0]++;
  193.  
  194.     if (result_mant[0] & 0x00200000)    /* overflow? */
  195.         {
  196.         num_exp--;
  197.         }
  198.     }
  199.   /* compute the resultant exponent */
  200.   num_exp = num_exp - den_exp + 0x3FF;
  201.  
  202.   /* reconstruct the result in local_num */
  203.   num_bdp->sign = result_sign;
  204.   num_bdp->exp = num_exp;
  205.   num_bdp->mant1 = result_mant[0] & 0xFFFFF;
  206.   num_bdp->mant2 = result_mant[1];
  207.   
  208.   /* done! */
  209.   return(local_num);
  210. }
  211. #endif
  212.  
  213. #ifdef L_muldf3
  214. /*double _muldf3 (a, b) double a, b; { return a * b; } */
  215.  
  216. /* a set of totally local kludges, for summing partial results into
  217.    result vector */
  218.  
  219. static unsigned long constant_zero_kludge = 0;
  220. #define XXADDL(partial, target) \
  221.     { register unsigned long * zero = &constant_zero_kludge + 1; \
  222.       asm volatile("addl  %3,%0@;\
  223.             addxl %1@-,%0@-" \
  224.             : "=a" (target), "=a" (zero) \
  225.             : "0" (target), "g" (partial), "1" (zero)); \
  226.     }
  227.  
  228. static char component_index[] = 
  229.     {
  230.     3, 3,                /* last ones */
  231.  
  232.     2, 3,                /* next to last x last */
  233.     3, 2,                /* ... */
  234.  
  235.     1, 3,
  236.     2, 2,
  237.     3, 1,
  238.  
  239.     0, 3,
  240.     1, 2,
  241.     2, 1,
  242.     3, 0,
  243.     
  244.     0, 2,
  245.     1, 1,
  246.     2, 0,
  247.     
  248.     0, 1,
  249.     1, 0,
  250.     
  251.     0, 0
  252.     };
  253.  
  254. qmult(m1_hi, m1_lo, m2_hi, m2_lo, result_hi, result_lo)
  255. unsigned long m1_hi, m1_lo, m2_hi, m2_lo, * result_hi, * result_lo;
  256. {
  257.   unsigned short * m1 = (unsigned short * )(&m1_hi);
  258.   unsigned short * m2 = (unsigned short * )(&m2_hi);
  259.   unsigned short result[10];        /* two extra for XADDL */
  260.   register unsigned long mult_register;
  261.   register unsigned long res1, res2, res3;
  262.   long * kludge;
  263.   short i, j;
  264.   char * idx_p = (char * )&component_index;
  265. /*
  266. fprintf(stderr, "  qmult: %08X:%08X * %08X:%08X\n", 
  267.     m1_hi, m1_lo, m2_hi, m2_lo);
  268. */
  269.   for (i = 0 ; i < 10 ; i++) result[i] = 0;
  270.  
  271. /* walk thru 'vectors' of mantissa pieces, doing unsigned multiplies
  272.    and summing results into result vector.  Note that the order is 
  273.    chosen to guarantee that no more adds than generated by XADDL are
  274.    needed.  */
  275.  
  276.   for ( ; ; )
  277.     {
  278.     i = *idx_p++;
  279.     j = *idx_p++;
  280.     mult_register = m1[i]; 
  281.     MUL(m2[j], mult_register);
  282.     kludge = (long * )(&result[i + j + 2]);
  283.     XXADDL(mult_register, kludge);
  284. /*
  285. fprintf(stderr, "  qmult: %d[%04X] * %d[%04X] -> %08X\n",
  286.   i, m1[i], j, m2[j], mult_register);
  287. fprintf(stderr, "       %04X %04X %04X %04X %04X %04X %04X %04X\n",
  288.   result[2], result[3], result[4], result[5], 
  289.   result[6], result[7], result[8], result[9]);
  290. */
  291.     if ((i == 0) && (j == 0))
  292.         break;
  293.     }
  294.  
  295.   /* get the result shifted to the right place.  Unfortunately, we need
  296.      the 53 bits that's 22 bits down in the result vec.  sigh */
  297.   kludge = (long * )(&result[2]);
  298.   res1 = *kludge++;
  299.   res2 = *kludge++;
  300.   res3 = *kludge;
  301.   for (i = 0 ; i < 12 ; i++)
  302.     ASL3(res1, res2, res3);
  303.   /* now put the result back */
  304.   *result_hi = res1;
  305.   *result_lo = res2;
  306. }
  307.  
  308. double dmult(f1, f2)
  309. double f1, f2;
  310. {
  311.     register long d2;
  312.     register unsigned d3, d4, d5, d6, d7;
  313.     unsigned long p1, p2;
  314.     
  315.     struct bitdouble
  316.     *bdp1 = (struct bitdouble *)&f1,
  317.     *bdp2 = (struct bitdouble *)&f2;
  318.     int exp1, exp2, i;
  319.     
  320.     exp1 = bdp1->exp; exp2 = bdp2->exp;
  321.     /* check for zero */
  322.     if (! exp1) return(0.0);
  323.     if (! exp2) return(0.0);
  324.     d2 = 0x00100000 | bdp1->mant1;
  325.     d3 = bdp1->mant2;
  326.     d4 = 0x00100000 | bdp2->mant1;
  327.     d5 = bdp2->mant2;
  328.     __qmult(d2, d3, d4, d5, &p1, &p2);
  329.     d6 = p1; d7 = p2;
  330.     
  331.     if (d6 & 0x00200000) {
  332.     ASR2(d6, d7);
  333.     exp1++;
  334.     }
  335.     
  336.     if (bdp1->sign == bdp2->sign) bdp1->sign = 0;
  337.     else bdp1->sign = 1;
  338.     bdp1->exp = exp1 + exp2 - 0x3ff;
  339.     bdp1->mant1 = d6;
  340.     bdp1->mant2 = d7;
  341.     return(f1);
  342. }
  343. #endif
  344.  
  345. #ifdef L_negdf2
  346. /*double _negdf2 (a) double a; { return -a; } */
  347.  
  348. double dneg(num)
  349. double num;
  350. {
  351.     long *i = (long *)#
  352.     if (*i)            /* only negate if non-zero */
  353.         *i ^= 0x80000000;
  354.     return(num);
  355. }
  356. #endif
  357.  
  358. #ifdef L_adddf3
  359. /*double _adddf3 (a, b) double a, b; { return a + b; } */
  360.  
  361. double dadd(f1, f2)
  362. double f1, f2;
  363. {
  364.     
  365.     register long d4, d5, d6, d7;
  366.     struct bitdouble
  367.     *bdp1 = (struct bitdouble *)&f1,
  368.         *bdp2 = (struct bitdouble *)&f2;
  369.     short exp1, exp2, sign1, sign2, howmuch, temp;
  370.     
  371.     exp1 = bdp1->exp; exp2 = bdp2->exp;
  372.     
  373.     /* check for zero */
  374.     if (! exp1) return(f2); if (! exp2) return(f1);
  375.     
  376.     /* align them */
  377.     if (exp1 < exp2) {
  378.     bdp1 = (struct bitdouble *)&f2; bdp2 = (struct bitdouble *)&f1;
  379.     exp1 = bdp1->exp; exp2 = bdp2->exp;
  380.     }
  381.     howmuch = exp1 - exp2;
  382.     if (howmuch > 53) return(f1);
  383.     
  384.     d7 = bdp2->mant1 | 0x00100000;
  385.     d6 = bdp2->mant2;
  386.     
  387.     d5 = bdp1->mant1 | 0x00100000;
  388.     d4 = bdp1->mant2;
  389.     
  390.     for (temp = 0; temp < howmuch; temp++) ASR2(d7, d6);
  391.     
  392.     /* take care of negative signs */
  393.     if (bdp1->sign) 
  394.     {
  395.     NEG(d5, d4);
  396.     }
  397.     if (bdp2->sign)
  398.     {
  399.     NEG(d7, d6);
  400.     }
  401.     
  402.     ADD2(d7, d6, d5, d4);
  403.     bdp1 = (struct bitdouble *)&f1;
  404.     
  405.     if (d5 < 0) {
  406.     NEG(d5, d4);
  407.     bdp1->sign = 1;
  408.     } else bdp1->sign = 0;
  409.     
  410.     if (d5 == 0 && d4 == 0) return(0.0);
  411.     
  412.     for(howmuch = 0; d5 >= 0; howmuch++) ASL2(d5, d4);
  413.     
  414.     ASL2(d5, d4);
  415.     for (temp = 0; temp < 12; temp++) ASR2(d5, d4);
  416.     bdp1->mant1 = d5;
  417.     bdp1->mant2 = d4;
  418.     bdp1->exp = exp1 + 11 - howmuch;
  419.     return(f1); 
  420. }
  421.  
  422. #endif
  423.  
  424. #ifdef L_subdf3
  425. double
  426. #ifdef __GCC_OLD__
  427.     _subdf3 (a, b)
  428. #else
  429.     __subdf3 (a, b)
  430. #endif
  431. double a, b;
  432. {
  433.     return a+(-b);
  434. }
  435. #endif
  436.  
  437. #ifdef L_cmpdf2
  438. /*
  439.   int _cmpdf2 (a, b) double a, b; { if (a > b) return 1;
  440.   else if (a < b) return -1; return 0; } 
  441.   */
  442.  
  443. int dcmp(f1, f2)
  444. double f1, f2;
  445. {
  446.     struct bitdouble *bdp1, *bdp2;
  447.     unsigned int s1, s2;
  448.     bdp1 = (struct bitdouble *)&f1;
  449.     bdp2 = (struct bitdouble *)&f2;
  450.     s1 = bdp1->sign;
  451.     s2 = bdp2->sign;
  452.     if (s1 > s2) return(-1);
  453.     if (s1 < s2) return(1);
  454.     /* if sign of both negative, switch them */
  455.     if (s1 != 0) {
  456.     bdp1 = (struct bitdouble *)&f2;
  457.     bdp2 = (struct bitdouble *)&f1;
  458.     }
  459.     s1 = bdp1->exp;
  460.     s2 = bdp2->exp;
  461.     if (s1 > s2) return(1);
  462.     if (s1 < s2) return(-1);
  463.     /* same exponent -- have to compare mantissa */
  464.     s1 = bdp1->mant1;
  465.     s2 = bdp2->mant1;
  466.     if (s1 > s2) return(1);
  467.     if (s1 < s2) return(-1);
  468.     s1 = bdp1->mant2;
  469.     s2 = bdp2->mant2;
  470.     if (s1 > s2) return(1);
  471.     if (s1 < s2) return(-1);
  472.     return(0); /* the same! */
  473. }
  474. #endif
  475.  
  476. #ifdef L_fixunsdfsi
  477. /*_fixunsdfsi (a) double a; { return (unsigned int) a; } */
  478.  
  479. unsigned long dtoul(f)
  480. double f;
  481. {
  482.     struct bitdouble *bdp;
  483.     int si, ex;
  484.     unsigned long mant1, mant2;
  485.     bdp = (struct bitdouble *)&f;
  486.     si = bdp->sign;
  487.     ex = bdp->exp;
  488.     mant1 = bdp->mant1 | 0x00100000;
  489.     mant2 = bdp->mant2;
  490.     
  491.     /* zero value */
  492. /*    if (ex == 0) return(0L); */ /* covered by next test */
  493.     /* less than 1 */
  494.     if (ex < 0x3ff) return(0L);
  495.     /* less than 0 */
  496.     if (si ) return(0L);
  497.     mant1 <<= 10;
  498.     mant1 += (mant2 >> 22);
  499.     mant1 >>= 30 + (0x3ff - ex);
  500.     return(mant1);
  501. }
  502.  
  503. long dtol(f)
  504. double f;
  505. {
  506.     struct bitdouble *bdp = (struct bitdouble *)&f;
  507.     
  508.     if (bdp->sign) {
  509.     bdp->sign = 0;
  510.     return( - dtoul(f));
  511.     }
  512.     return(dtoul(f));
  513. }
  514. #endif
  515.  
  516. #ifdef L_fixunsdfdi
  517.  
  518. /*double _fixunsdfdi (a) double a; { union double_di u;
  519.   u.i[LOW] = (unsigned int) a; u.i[HIGH] = 0; return u.d; } */
  520. #endif
  521.  
  522. #ifdef L_fixdfdi
  523. double
  524. #ifdef __GCC_OLD__
  525.     _fixdfdi (a)
  526. #else
  527.     __fixdfdi (a)
  528. #endif
  529. double a;
  530. {
  531.     union double_di u;
  532.     u.i[LOW] = (int) a;
  533.     u.i[HIGH] = (int) a < 0 ? -1 : 0;
  534.     return u.d;
  535. }
  536. #endif
  537.  
  538. #ifdef L_floatsidf
  539. /* double _floatsidf (a) int a; { return (double) a; } */
  540.  
  541. double ltod(i)
  542. long i;
  543. {
  544.     int expo, shift;
  545.     double retval;
  546.     struct bitdouble *bdp = (struct bitdouble *)&retval;
  547.     if (i == 0) {
  548.     long *t = (long *)&retval;
  549.     t[0] = 0L;
  550.     t[1] = 0L;
  551.     return(retval);
  552.     }
  553.     if (i < 0) {
  554.     bdp->sign = 1;
  555.     i = -i;
  556.     } else bdp->sign = 0;
  557.     shift = i;
  558.     for (expo = 0x3ff + 31 ; shift > 0; expo--, shift <<= 1);
  559.     shift <<= 1;
  560.     bdp->exp = expo;
  561.     bdp->mant1 = shift >> 12;
  562.     bdp->mant2 = shift << 20;
  563.     return(retval);
  564. }
  565.  
  566. #endif
  567.  
  568. #ifdef L_floatdidf
  569. /* ok as is -- will call other routines */
  570. double
  571. #ifdef __GCC_OLD__
  572.     _floatdidf (u)
  573. #else
  574.     __floatdidf (u)
  575. #endif
  576. union double_di u;
  577. {
  578.     register double hi
  579.     = ((double) u.i[HIGH]) * (double) 0x10000 * (double) 0x10000;
  580.     register double low = (unsigned int) u.i[LOW];
  581.     return hi + low;
  582. }
  583. #endif
  584.  
  585. #ifdef L_truncdfsf2
  586. float __dtof(d)
  587. double d;
  588. {
  589.     struct bitdouble *bdp = (struct bitdouble *)&d;
  590.     float retval;
  591.     struct bitfloat *bfp = (struct bitfloat *)&retval;
  592.     int tempval;
  593.     
  594.     bfp->sign = bdp->sign;
  595.     if (bdp->exp == 0) return ((float) 0.0);
  596.     bfp->exp = bdp->exp - 0x400 + 0x80;
  597.     tempval = (bdp->mant1 << 4 ) + ((0xF0000000 & bdp->mant2) >> 28);
  598.     /* round */
  599.     tempval++;
  600.     if (tempval == 0x01000000) bfp->exp++;
  601.     bfp->mant = tempval >> 1;
  602.     return(retval);
  603. }
  604.  
  605. SFVALUE
  606. #ifdef __GCC_OLD__
  607.     _truncdfsf2 (a)
  608. #else
  609.     __truncdfsf2 (a)
  610. #endif
  611. double a;
  612. {
  613.     union flt_or_int intify;
  614.     return INTIFY (__dtof(a));
  615. }
  616. #endif
  617.  
  618. #ifdef L_extendsfdf2
  619. double __ftod(f)
  620. union flt_or_int f;
  621. {
  622.     double retval;
  623.     struct bitfloat *bfp = (struct bitfloat *)&f.f;
  624.     struct bitdouble *bdp = (struct bitdouble *)&retval;
  625.     if (bfp->exp == 0) return(0.0);
  626.     bdp->sign = bfp->sign;
  627.     bdp->exp = 0x400 - 0x80 + bfp->exp;
  628.     bdp->mant1 = bfp->mant >> 3;
  629.     bdp->mant2 = (bfp->mant & 0x7) << 29;
  630.     /*printf("returning %f from extendsfdf2\n", retval);*/
  631.     return(retval);
  632. }
  633.  
  634. double
  635. #ifdef __GCC_OLD__
  636.     _extendsfdf2 (a)
  637. #else
  638.     __extendsfdf2 (a)
  639. #endif
  640. union flt_or_int a;
  641. {
  642.     union flt_or_int intify;
  643.     double retval;
  644.     return (__ftod(a));
  645. }
  646. #endif
  647.